home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / lib / mathlib / libconv / TRY / check1d.c next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  15.6 KB  |  719 lines
  1. #include <stdio.h>
  2. #include <sys/time.h>
  3. #include <math.h>
  4. #include "conv.h"
  5.  
  6. /* --------- The following definitions change 
  7.         according to precision required -------- */
  8.  
  9.  
  10. #ifdef SINGLE       /* complex single precision    */
  11.  
  12. typedef float  this_variable;
  13. typedef float  this_real;
  14. #define MAX_RECUR    17
  15. #define OPS_PER_ITER    2
  16.  
  17. #define SIMPLE_FIR        simple_sfir1d_
  18. #define ORNL_FIR        ornl_sfir1d_
  19. #define MY_FIR            sfir1d_
  20. #define MY_C_FIR        sfir1d
  21.  
  22. #define SIMPLE_IIR        simple_siir1d_
  23. #define ORNL_IIR        ornl_siir1d_
  24. #define MY_IIR            siir1d_
  25. #define MY_C_IIR        siir1d
  26.  
  27. #define SIMPLE_COR        simple_scor1d_
  28. #define MY_COR            scor1d_
  29. #define MY_INIT            sinit_
  30. #define TOLERANCE        1.e-3
  31. #define THIS_IS_REAL
  32. #endif
  33. #ifdef DOUBLE            /* real double precision    */
  34.  
  35. typedef double  this_variable;
  36. typedef double  this_real;
  37.  
  38. #define MAX_RECUR    17
  39. #define OPS_PER_ITER    2
  40.  
  41. #define SIMPLE_FIR        simple_dfir1d_
  42. #define ORNL_FIR        ornl_dfir1d_
  43. #define MY_FIR            dfir1d_
  44. #define MY_C_FIR        dfir1d
  45.  
  46. #define SIMPLE_IIR        simple_diir1d_
  47. #define ORNL_IIR        ornl_diir1d_
  48. #define MY_IIR            diir1d_
  49. #define MY_C_IIR        diir1d
  50.  
  51. #define SIMPLE_COR        simple_dcor1d_
  52. #define MY_COR            dcor1d_
  53. #define MY_INIT            dinit_
  54. #define TOLERANCE        1.e-7
  55. #define THIS_IS_REAL
  56.  
  57. #endif
  58. #ifdef COMPLEX            /* complex single precision    */
  59.  
  60. typedef complex this_variable;
  61. typedef float  this_real;
  62.  
  63. #define MAX_RECUR    17
  64. #define OPS_PER_ITER    8
  65.  
  66. #define SIMPLE_FIR        simple_cfir1d_
  67. #define ORNL_FIR        ornl_cfir1d_
  68. #define MY_FIR            cfir1d_
  69. #define MY_C_FIR        cfir1d
  70.  
  71. #define SIMPLE_IIR        simple_ciir1d_
  72. #define ORNL_IIR        ornl_ciir1d_
  73. #define MY_IIR            ciir1d_
  74. #define MY_C_IIR        ciir1d
  75.  
  76. #define SIMPLE_COR        simple_ccor1d_
  77. #define MY_COR            ccor1d_
  78. #define MY_INIT            cinit_
  79. #define TOLERANCE        1.e-3
  80. #define THIS_IS_COMPLEX
  81.  
  82. #endif
  83. #ifdef ZOMPLEX            /* complex double precision    */
  84.  
  85. typedef zomplex this_variable;
  86. typedef double  this_real;
  87.  
  88. #define MAX_RECUR    17
  89. #define OPS_PER_ITER    8
  90.  
  91. #define SIMPLE_FIR        simple_zfir1d_
  92. #define ORNL_FIR        ornl_zfir1d_
  93. #define MY_FIR            zfir1d_
  94. #define MY_C_FIR        zfir1d
  95.  
  96. #define SIMPLE_IIR        simple_ziir1d_
  97. #define ORNL_IIR        ornl_ziir1d_
  98. #define MY_IIR            ziir1d_
  99. #define MY_C_IIR        ziir1d
  100.  
  101. #define SIMPLE_COR        simple_zcor1d_
  102. #define MY_COR            zcor1d_
  103. #define MY_INIT            zinit_
  104. #define TOLERANCE        1.e-7
  105. #define THIS_IS_COMPLEX
  106.  
  107. #endif
  108.  
  109. /* ---------- The rest is the same    ---------------- */
  110.  
  111. void MY_FIR(), SIMPLE_FIR(), ORNL_FIR();
  112. void MY_IIR(), MY_C_IIR(), SIMPLE_IIR(), ORNL_IIR();
  113.  
  114. #define MAX_VAL 3
  115. this_real my_val[MAX_VAL] = {0., -.33, 1.};
  116.  
  117. #define MAX_SIZE        35
  118. #define DEFAULT_TRIALS        999
  119. #define MAX_INC            3
  120. #define MAX_TIMES        3
  121.  
  122.  
  123. #define ABS(a) ( ((a)>0) ? (a) : -(a))
  124. #define MAX(a,b) ( ((a) > (b)) ? (a) : (b))
  125. #define MIN(a,b) ( ((a) < (b)) ? (a) : (b))
  126.  
  127. void (*simple_function)();
  128. void (*ornl_function)();
  129. void (*my_function)();
  130. void GetArguments();
  131.  
  132. double check_this();
  133.  
  134. int parallel;
  135. int is_random;
  136. int n_trials;
  137. int len = 4;
  138.  
  139.     int size, n_trials, n_times, nf, ng, nh;
  140.     int incf, nf1, nf2, incg, ng1, ng2, inch, nh1, nh2;
  141.     this_variable *vx, *vy, *vz, *vt, *vr;
  142.     this_variable alpha, beta, one, zero;
  143.     double t, x, y, z, u;
  144.  
  145. main(argc,argv)
  146. int argc;
  147. char *argv[];
  148. {
  149.     int i, j, k;
  150.  
  151. /* ******************************************************* */
  152.     GetArguments( argc, argv);
  153. /* ******************************************************* */
  154.  
  155.   for( ; n_trials > 0 ; n_trials --) { 
  156.  
  157.     get_values();
  158.  
  159.     printf("%3dx( %2d,%3d,%3d | %2d,%3d,%3d | %2d,%3d,%3d <>",
  160.         n_times, incf,nf1,nf, incg,ng1,ng, inch,nh1,nh);
  161. #ifdef THIS_IS_REAL
  162.     printf(" %5.2f, %5.2f ):", alpha, beta);
  163. #endif
  164. #ifdef THIS_IS_COMPLEX
  165.     printf(" [%5.2f,%5.2f], [%5.2f,%5.2f] ):", 
  166.         alpha.real, alpha.imag, beta.real, beta.imag);
  167. #endif
  168.  
  169.     fflush(stdout);
  170.  
  171.         do_it();
  172.  
  173.     printf("\n", x);
  174.       }
  175.     return(0);
  176. }
  177. do_it() 
  178. {
  179.     int i,j,ii, jj, k, max_size;
  180.     int this_inc, that_inc, i1, i2, j1, j2, k1, k2, l1, l2;
  181.     this_variable *this_vx, *this_vy, *this_vz, *this_vt, *this_vr;
  182.  
  183.     vx = (this_variable *)my_malloc(((nf+1)*MAX_INC) * sizeof( this_variable));
  184.     vy = (this_variable *)my_malloc(((ng+1)*MAX_INC) * sizeof( this_variable));
  185.     vr = (this_variable *)my_malloc(((ng+1)*MAX_INC) * sizeof( this_variable));
  186.     vz = (this_variable *)my_malloc(((nh+1)*MAX_INC) * sizeof( this_variable));
  187.     vt = (this_variable *)my_malloc(((nh+1)*MAX_INC) * sizeof( this_variable));
  188.  
  189.     if( (vx == (this_variable *)0) || 
  190.     (vy == (this_variable *)0) ||
  191.     (vr == (this_variable *)0) ||
  192.     (vt == (this_variable *)0) ||
  193.     (vz == (this_variable *)0))  {
  194.     fprintf( stderr, "Malloc problem ... Exiting");
  195.     exit( 0 );
  196.     }
  197.     i = (nf+1)*MAX_INC;
  198.     i = MY_INIT( &i, vx);
  199.     i = (ng+1)*MAX_INC;
  200.     j = MY_INIT( &i, vy);
  201.  
  202.     max_size = (nh+1)*MAX_INC;
  203. /*
  204. *   Checking against a simple FIR convolution
  205. */
  206.     MY_INIT( &max_size, vz);
  207.     bcopy( vz, vt, max_size*sizeof(this_variable));
  208.  
  209.     for( ii = 0; ii < n_times ; ii ++)
  210.         SIMPLE_FIR( &nf, &ng, vx, vy, vz);
  211.  
  212.     this_inc = 1;
  213.     i1 = 0; i2 = nf;
  214.     j1 = 0; j2 = ng;
  215.     k1 = 0; k2 = i2+j2-1;
  216.     for( ii = 0; ii < n_times ; ii ++)
  217.         MY_FIR(    vx,&this_inc,&i1,&i2,
  218.             vy,&this_inc,&j1,&j2,
  219.             vt,&this_inc,&k1,&k2,
  220.             &one, &zero);
  221.  
  222.     t = check_this( vz, vt, max_size) ;
  223.  
  224. /*
  225. *   Checking against a in-place convolution with simple FIR convolution
  226. */
  227.     MY_INIT( &max_size, vz);
  228.     bcopy( vz, vt, max_size*sizeof(this_variable));
  229.     bcopy( vx, vt, (nf)*sizeof(this_variable));
  230.  
  231.         SIMPLE_FIR( &nf, &ng, vx, vy, vz);
  232.  
  233.     this_inc = 1;
  234.     i1 = 0; i2 = nf;
  235.     j1 = 0; j2 = ng;
  236.     k1 = 0; k2 = i2+j2-1;
  237.  
  238.         MY_FIR(    vt,&this_inc,&i1,&i2,
  239.             vy,&this_inc,&j1,&j2,
  240.             vt,&this_inc,&k1,&k2,
  241.             &one, &zero);
  242.  
  243.     t = check_this( vz, vt, max_size) ;
  244.  
  245. /*
  246. *   Checking against a simple IIR convolution
  247. */
  248.     MY_INIT( &max_size, vz);
  249.     bcopy( vz, vt, max_size*sizeof(this_variable));
  250.  
  251.     this_inc = 1;
  252.     i1 = 0; i2 = nf;
  253.     j1 = 0; j2 = ng;
  254.     k1 = 0; k2 = nf;
  255.     init_min_phase( j2, vr, this_inc);
  256.  
  257.     for( ii = 0; ii < n_times ; ii ++)
  258.         SIMPLE_IIR( &nf, &j2, vx, vr, vz);
  259.  
  260.     for( ii = 0; ii < n_times ; ii ++)
  261.         MY_IIR(        vx,&this_inc,&i1,&i2,
  262.             vr,&this_inc,&j1,&j2,
  263.             vt,&this_inc,&k1,&k2);
  264.  
  265.     t = check_this( vz, vt, max_size) ;
  266. /*
  267. *   Checking against a in-place convolution with simple IIR convolution
  268. */
  269.     max_size = (nh+1)*MAX_INC;
  270.     MY_INIT( &max_size, vz);
  271.     bcopy( vz, vt, (max_size)*sizeof(this_variable));
  272.     bcopy( vx, vt, (nf)*sizeof(this_variable));
  273.  
  274.     this_inc = 1;
  275.     i1 = 0; i2 = nf;
  276.     j1 = 0; j2 = ng;
  277.     k1 = 0; k2 = i2;
  278.  
  279.     init_min_phase( j2, vr, this_inc);
  280.  
  281.     SIMPLE_IIR( &i2, &j2, vx, vr, vz);
  282.  
  283.     MY_IIR(        vt,&this_inc,&i1,&i2,
  284.             vr,&this_inc,&j1,&j2,
  285.             vt,&this_inc,&k1,&k2);
  286.  
  287.     t = check_this( vz, vt, max_size) ;
  288. /*
  289. *   Checking FIR against IIR  ...  IIR*IIR = ident
  290. */
  291.     MY_INIT( &max_size, vz);
  292.     bcopy( vz, vt, max_size*sizeof(this_variable));
  293.  
  294.     this_inc = 1;
  295.     i1 = 0; i2 = MIN(nf, MAX_RECUR);
  296.     j1 = 0; j2 = MIN(ng, MAX_RECUR);
  297.     k1 = 0; k2 = MIN(nf, MAX_RECUR);
  298.  
  299.     init_min_phase( j2, vr, this_inc);
  300.  
  301.     MY_FIR(        vt,&this_inc,&i1,&i2,
  302.             vr,&this_inc,&j1,&j2,
  303.             vt,&this_inc,&k1,&k2,
  304.             &one, &zero);
  305.  
  306.     MY_IIR(        vt,&this_inc,&k1,&k2,
  307.             vr,&this_inc,&j1,&j2,
  308.             vt,&this_inc,&k1,&k2);
  309.  
  310.     t = check_this( vz, vt, max_size) ;
  311. /*
  312. *   Checking COR against a simple correlation
  313. */
  314.     bcopy( vz, vt, (max_size)*sizeof(this_variable));
  315.  
  316.     for( ii = 0; ii < n_times ; ii ++)
  317.     SIMPLE_COR( &nf, &ng, vx, vy, vz);
  318.     
  319.     this_inc = 1;
  320.     i1 = 0; i2 = nf;
  321.     j1 = 0; j2 = ng;
  322.     k1 = i1-(j2-1); k2 = nf+ng-1;
  323.     for( ii = 0; ii < n_times ; ii ++)
  324.     MY_COR(        vx,&this_inc,&i1,&i2,
  325.             vy,&this_inc,&j1,&j2,
  326.             vt,&this_inc,&k1,&k2);
  327.  
  328.     t = check_this( vz, vt, max_size) ;
  329. /*
  330. *   Checking FIR against a simple correlation
  331. */
  332.     bcopy( vz, vt, (max_size)*sizeof(this_variable));
  333.  
  334.     for( ii = 0; ii < n_times ; ii ++)
  335.     SIMPLE_COR( &nf, &ng, vx, vy, vz);
  336.     
  337.     this_inc = 1;
  338.     i1 = 0; i2 = nf;
  339.     j1 = -(ng-1); j2 = ng;
  340.     k1 = i1+j1; k2 = i2+j2-1;
  341.     this_vy = vy + j2-1;
  342.     that_inc = -1;
  343.     for( ii = 0; ii < n_times ; ii ++)
  344.     MY_FIR(        vx,&this_inc,&i1,&i2,
  345.             this_vy,&that_inc,&j1,&j2,
  346.             vt,&this_inc,&k1,&k2,
  347.             &one, &zero);
  348.  
  349.     t = check_this( vz, vt, max_size) ;
  350. /*
  351. *   Checking FIR against a Fortran template
  352. */
  353.     this_vx = ( incf < 0) ? (vx - incf * (nf - 1)) : vx;
  354.     this_vy = ( incg < 0) ? (vy - incg * (ng - 1)) : vy;
  355.     this_vr = ( incg < 0) ? (vr - incg * (ng - 1)) : vr;
  356.     this_vz = ( inch < 0) ? (vz - inch * (nh - 1)) : vz;
  357.     this_vt = ( inch < 0) ? (vt - inch * (nh - 1)) : vt;
  358.  
  359.     i1 = nf1; i2 = nf;
  360.     j1 = ng1; j2 = ng;
  361.     k1 = nh1; k2 = nh;
  362.  
  363.     MY_INIT( &max_size, vz);
  364.     bcopy( vz, vt, max_size * sizeof(this_variable));
  365.  
  366.     for( ii = 0; ii < n_times ; ii ++)
  367.     ORNL_FIR(   this_vx,&incf,&i1,&i2,
  368.             this_vy,&incg,&j1,&j2,
  369.             this_vz,&inch,&k1,&k2,
  370.             &alpha, &beta);
  371.  
  372.     for( ii = 0; ii < n_times ; ii ++)
  373.     MY_FIR(        this_vx,&incf,&i1,&i2,
  374.             this_vy,&incg,&j1,&j2,
  375.             this_vt,&inch,&k1,&k2,
  376.             &alpha, &beta);
  377.  
  378.     t = check_this( vz, vt, max_size) ;
  379. /*
  380. *   Checking the C_FIR version against the Fortran template
  381. */
  382.     MY_INIT( &max_size, vz);
  383.     bcopy( vz, vt, (max_size)*sizeof(this_variable));
  384.  
  385.     for( ii = 0; ii < n_times ; ii ++)
  386.     ORNL_FIR(   this_vx,&incf,&i1,&i2,
  387.             this_vy,&incg,&j1,&j2,
  388.             this_vz,&inch,&k1,&k2,
  389.             &alpha, &beta);
  390.  
  391.     for( ii = 0; ii < n_times ; ii ++)
  392.     MY_C_FIR(   this_vx,incf,i1,i2,
  393.             this_vy,incg,j1,j2,
  394.             this_vt,inch,k1,k2,
  395. #ifdef THIS_IS_REAL
  396.             alpha, beta);
  397. #endif
  398. #ifdef THIS_IS_COMPLEX
  399.             &alpha, &beta);
  400. #endif
  401.  
  402.     t = check_this( vz, vt, max_size) ;
  403. /*
  404. *   Checking the IIR against the Fortran template
  405. */
  406.     MY_INIT( &max_size, vz);
  407.     bcopy( vz, vt, (max_size)*sizeof(this_variable));
  408.  
  409.     init_min_phase( j2, this_vr, incg);
  410.  
  411.     ORNL_IIR(   this_vx,&incf,&i1,&i2,
  412.             this_vr,&incg,&j1,&j2,
  413.             this_vt,&inch,&k1,&k2);
  414.  
  415.     MY_IIR(        this_vx,&incf,&i1,&i2,
  416.             this_vr,&incg,&j1,&j2,
  417.             this_vz,&inch,&k1,&k2);
  418.  
  419.     t = check_this( vz, vt, max_size) ;
  420. /*
  421. *   Checking the C_IIR against the Fortran template
  422. */
  423.     MY_INIT( &max_size, vz);
  424.     bcopy( vz, vt, (max_size)*sizeof(this_variable));
  425.  
  426.     ORNL_IIR(   this_vx,&incf,&i1,&i2,
  427.             this_vr,&incg,&j1,&j2,
  428.             this_vt,&inch,&k1,&k2);
  429.  
  430.     MY_C_IIR(   this_vx,incf,i1,i2,
  431.             this_vr,incg,j1,j2,
  432.             this_vz,inch,k1,k2);
  433.  
  434.     t = check_this( vz, vt, max_size) ;
  435. /*
  436. *
  437. */
  438.     my_free ( vt);
  439.     my_free ( vz);
  440.     my_free ( vr);
  441.     my_free ( vy);
  442.     my_free ( vx);
  443. }
  444.  
  445. void GetArguments( argc, argv)
  446. int argc;
  447. char *argv[];
  448. {
  449.     int i, j, k;
  450.     int nerror = 0;
  451.  
  452. #define ON    1
  453.  
  454.     parallel = 0;
  455.     is_random = 1;
  456.     n_trials = DEFAULT_TRIALS;
  457.  
  458. /* ******************************************************* */
  459.     for ( i = 1 ; (i < argc) && (nerror != ON) ; i ++ ) {
  460.     if( argv[i][0] == '-') {
  461.         switch ( argv[i][1]) {
  462.  
  463.         case 'i' :
  464.         case 'I' :  
  465.                 is_random = 0;
  466.                 break;
  467.         default  :  nerror = ON;
  468.         }
  469.     }
  470.     else {
  471.         if( i+1 > argc)
  472.         nerror = ON;
  473.         else { 
  474.         n_trials = atoi( argv[i]);
  475.         }
  476.     }
  477.     }
  478.     if( nerror == ON) {
  479.     fprintf( stderr, 
  480.         "Usage : %s [-p(arallel)] [n_trials]\n", argv[0]);
  481.     exit(0);
  482.     }
  483.     simple_function = SIMPLE_FIR;
  484.     ornl_function = ORNL_FIR;
  485.     my_function = MY_FIR;
  486.  
  487.     if( is_random)
  488.     srandom( (123*getpid()) | 0x01);
  489.     else 
  490.     srandom( 1);
  491. }
  492. get_values(){ 
  493.   int jj;
  494.     if( is_random) { 
  495.     nf = random() % MAX_SIZE + 1;
  496.     ng = random() % MAX_SIZE + 1;
  497.     nh = nf+ng+random() % MAX_SIZE;
  498.  
  499.     nf1 = random() % nf;
  500.     nf2 = nf1 + random() % (nf - nf1);
  501.     ng1 = random() % ng;
  502.     ng2 = ng1 + random() % (ng - ng1);
  503.     nh1 = nf1+ng1;
  504.     nh2 = nh1 + random() % (nh - nh1);
  505.  
  506.     jj = random() % 23 - 11;
  507.     nf1 -= jj; nf2 -= jj;
  508.     jj = random() % 23 - 11;
  509.     ng1 -= jj; ng2 -= jj;
  510.     jj = random() % 23 - 11;
  511.     nh1 += jj; nh2 += jj;
  512.     
  513.     incf = random() % (2 * MAX_INC) - MAX_INC;
  514.     incg = random() % (2 * MAX_INC) - MAX_INC;
  515.     inch = random() % (2 * MAX_INC) - MAX_INC;
  516.  
  517.     if(incf == 0) incf = 1;
  518.     if(incg == 0) incg = 1;
  519.     if(inch == 0) inch = 1;
  520.  
  521.     n_times = random() % MAX_TIMES + 1;
  522. #ifdef THIS_IS_REAL
  523.     alpha = my_val[ random() % MAX_VAL];
  524.     beta  = my_val[ random() % MAX_VAL];
  525. #endif
  526. #ifdef THIS_IS_COMPLEX
  527.     alpha.real = my_val[ random() % MAX_VAL];
  528.     alpha.imag = my_val[ random() % MAX_VAL];
  529.     beta.real = my_val[ random() % MAX_VAL];
  530.     beta.imag = my_val[ random() % MAX_VAL];
  531. #endif
  532.     }
  533.     else {
  534.     printf( "Enter sizex : ");
  535.     scanf( "%d", &nf);
  536.     printf( "Enter sizey : ");
  537.     scanf( "%d", &ng);
  538.     printf( "Enter sizez : ");
  539.     scanf( "%d", &nh);
  540.  
  541.     n_times = 1;
  542.  
  543.     printf( "Enter x0 : ");
  544.     scanf( "%d", &nf1);
  545.     nf2 = nf1+nf-1;
  546.  
  547.     printf( "Enter y0 : ");
  548.     scanf( "%d", &ng1);
  549.     ng2 = ng1+ng-1;
  550.  
  551.     printf( "Enter z0 : ");
  552.     scanf( "%d", &nh1);
  553.     nh2 = nh1+nh-1;
  554.  
  555.     printf( "Enter incx : ");
  556.     scanf( "%d", &incf);
  557.     printf( "Enter incy : ");
  558.     scanf( "%d", &incg);
  559.     printf( "Enter incz : ");
  560.     scanf( "%d", &inch);
  561.  
  562. #ifdef THIS_IS_REAL
  563.     alpha = .23;
  564.     beta  = -.41;
  565. #endif
  566. #ifdef THIS_IS_COMPLEX
  567.     alpha.real = .22;
  568.     alpha.imag = -.11;
  569.     beta.real = -.22;
  570.     beta.imag = .33;
  571. #endif
  572.     }
  573. #ifdef THIS_IS_REAL
  574.     one = 1.;
  575.     zero = 0.;
  576. #endif
  577. #ifdef THIS_IS_COMPLEX
  578.     one.real = 1.;
  579.     one.imag = 0.;
  580.     zero.real = 0.;
  581.     zero.imag = 0.;
  582. #endif
  583.  
  584. }
  585.  
  586. double check_this( this_variable *a, this_variable *b, int n)
  587. {
  588.     double max, av;
  589.     int i, nerr=0;
  590.     max = av = 0.0;
  591.  
  592. #ifdef THIS_IS_REAL
  593.     for (i = 0; i < n; i++){
  594.     t = a[i] - b[i];
  595.     t = t * t;
  596.     if( t > max) max = t;
  597.     t = a[i];
  598.     av += t * t;
  599.     }
  600. #endif
  601. #ifdef THIS_IS_COMPLEX
  602.     for (i = 0; i < n; i++){
  603.     t = a[i].real - b[i].real;
  604.     z = a[i].imag - b[i].imag;
  605.     t = t * t + z * z;
  606.     if( t > max) max = t;
  607.     t = a[i].real;
  608.     z = a[i].imag;
  609.     av += t * t + z * z;
  610.     }
  611. #endif
  612.     max = sqrt(max);
  613.     av  = sqrt( av / (double)n);
  614.     if( av > 0.0) {
  615.     max = max / av;
  616.     }
  617.     if( max > TOLERANCE) {
  618.     fprintf(stderr, " Difference with reference is %.16e \n\n", max);
  619.     max = TOLERANCE * av;
  620.     max = max * max;
  621.     for (i = 0; i < n; i++){
  622. #ifdef THIS_IS_REAL
  623.         t = a[i] - b[i];
  624.         t = t * t;
  625. #endif
  626. #ifdef THIS_IS_COMPLEX
  627.         t = a[i].real - b[i].real;
  628.         z = a[i].imag - b[i].imag;
  629.         t = t * t + z * z;
  630. #endif
  631.         if( t > max) {
  632.         if( nerr < 10)
  633.             fprintf(stderr, " %d", i);
  634.         nerr++;
  635.         }
  636.     }
  637.     fprintf( stderr, "\n %d Errors detected \n", nerr);
  638.  
  639.     exit(0);
  640.     }
  641.     else
  642.     printf("-", t);
  643.  
  644.      fflush(stdout);
  645.  
  646.     return(max);
  647. }
  648. int init_min_phase( int ng, this_variable *vr, int inc)
  649. {
  650. #define MAX_ROOT .4
  651.  
  652.     this_variable duplet[2], *this, *that, *tmp;
  653.     int i, ione=1, two = 2;
  654.  
  655.     this = (this_variable *) my_malloc( ng * sizeof(this_variable));
  656.     that = (this_variable *) my_malloc( ng * sizeof(this_variable));
  657.  
  658. #ifdef THIS_IS_REAL
  659.         *duplet = 1.0;
  660.         *this  = 1.0;
  661. #endif
  662. #ifdef THIS_IS_COMPLEX
  663.     (*duplet).real = 1.0;
  664.     (*duplet).imag = 0.0;
  665.     (*this).real = (*duplet).real;
  666.     (*this).imag = (*duplet).imag;
  667. #endif
  668.  
  669.     for (i = 1; i < ng ; i++) {
  670.     MY_INIT( &ione, duplet+1);
  671. /*
  672.     printf(" %d %f \n", i, duplet[1]);
  673. */
  674. #ifdef THIS_IS_REAL
  675.     duplet[1] *= MAX_ROOT;
  676. #endif
  677. #ifdef THIS_IS_COMPLEX
  678.     duplet[1].real *= MAX_ROOT;
  679.     duplet[1].imag *= MAX_ROOT;
  680. #endif
  681.     SIMPLE_FIR( &i, &two, this, duplet, that);
  682.  
  683.     tmp  = this;
  684.     this = that;
  685.     that = tmp;
  686.     }
  687.     for ( i = 0, tmp = vr; i < ng; i++, tmp += inc) {
  688.     *tmp = this[i];
  689.     }
  690.     my_free( that);
  691.     my_free (this);
  692. /*
  693. */
  694. #ifdef THIS_IS_REAL
  695.     if( vr[0] != 1.) { 
  696. #endif
  697. #ifdef THIS_IS_COMPLEX
  698.     if( (vr->real != 1.) || (vr->imag != 0.) ) {
  699. #endif
  700.     fprintf(stderr, "Error in First sample of Min_Phase filter \n");
  701.     }
  702.  
  703. #ifdef THIS_IS_REAL
  704.     if( ABS(vr[(ng-1)*inc]) > 1.) { 
  705.     fprintf(stderr, "Error in Last sample of Min_Phase filter: %.16e \n", vr[ng-1]);
  706.     }
  707. #endif
  708. #ifdef THIS_IS_COMPLEX
  709.     if( (ABS(vr[(ng-1)*inc].real) > 1.) || (ABS(vr[(ng-1)*inc].imag) > 1.) ) {
  710.     fprintf(stderr, "Error in Last sample of Min_Phase filter: (%.16e, %.16e) \n", 
  711.         vr[(ng-1)*inc].real, vr[(ng-1)*inc].imag);
  712.     }
  713. #endif
  714. #ifdef THIS_IS_COMPLEX
  715.     (*vr).imag = 0.5;
  716. #endif
  717.  
  718. }
  719.